Stacking Interactions in Denaturation of DNA Fragments 
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A mesoscopic model for heterogeneous DNA denaturation is developed in the framework of the 
path integral formalism. The base pair stretchings are treated as one-dimensional, time dependent 
paths contributing to the partition function. The size of the paths ensemble, which measures the 
degree of cooperativity of the system, is computed versus temperature consistently with the model 
potential physical requirements. It is shown that the ensemble size strongly varies with the molecule 
backbone stiffness providing a quantitative relation between stacking and features of the melting 
transition. The latter is an overall smooth crossover which begins from the adenine-thymine rich 
portions of the fragment. The harmonic stacking coupling shifts, along the T-axis, the occurrence 
of the multistep denaturation but it does not change the character of the crossover. The methods 
to compute the fractions of open base pairs versus temperature are discussed: by averaging the base 
pair displacements over the path ensemble we find that such fractions signal the multisteps of the 
transition in good agreement with the indications provided by the specific heat plots. 

PACS numbers: 87.14.gk, 87.15.A-, 87.15.Zg, 05.10.-a 



1. Introduction 



It has been long recognized that the physical properties 
of the DNA molecule are key to understand its biological 
function [l| . The core of the double helix embodies the 
genetic code through the sequence of base pairs. Gene 
transcription occurs as the hydrogen bonds between com- 
plementary bases can be broken thus leading to the for- 
mation of a transcription bubble in which the bases are 
exposed for chemical reactions Q. While in heteroge- 
neous fragments the bubble size is larger when the con- 
tent of weakly bonded adenine-thymine (AT) base pairs 
is higher, bubble formation can be also achieved exper- 
imentally by gradually heating the system as far as the 
two strands eventually separate. This process, known as 
thermal denaturation or melting, has been extensively 
studied in the last decades due to its relevance for the 
comprehension of the replication and transcription mech- 
anisms Q. However, agreement has not been reached 
so far regarding the character of the melting transition 
whether first- or second- order. Besides being intriguing 
for the statistical physicists community the question has 
biological importance also in view of recent investigations 
pointing to correlations between thermodynamical melt- 
ing properties and coding sequences in genomes 
Thermal denaturation is exploited in molecular biology 
where an understanding of the sequence specificity of 
melting is desirable for polymerase chain reaction. 

Two families of theoretical approaches have been devel- 
oped to account for the melting phenomenon. The first 
is based on the Poland-Scheraga model [8l4llj which sees 
the DNA as a two state Ising-like sequence of base pairs 
with regions of variable size, denaturated loops, opening 
temporarily due to thermal fluctuations. In this picture 
self-avoiding interactions between loops and the rest of 
the helix [12| may sharpen the denaturation leading to a 



first order transition in two dimensions and above [l3l Il4j 
while a mechanically induced denaturation should remain 
second order [l5j . Improvements of Ising-like models, in- 
troducing rotational degrees of freedom and accounting 
for bending rigidity of bubbles, reduce the importance of 
loop entropy and also suggest that the denaturation is 
rather a smooth crossover [161 ] . The second approach as- 
sumes the Peyrard-Bishop (PB) model [13] in which the 
stretching between complementary base pairs is repre- 
sented by a one dimensional, continuous variable bring- 
ing the advantage that also intermediate states in the 
DNA dynamics can be described. As a key refinement of 
the PB Hamiltonian, anharmonic stiffness has been in- 
troduced in the intra-strand stacking potential [l8[ thus 
accounting for an entropy driven sharp denaturation [l9| . 
While the PB model has proved successful to predict 
multistep melting in heterogeneous DNA [2(| [5l|, sev- 
eral methodologies such as molecular dynamics, transfer 
integral calculations [22j and renormalization group tech- 
niques [23j have been employed to investigate the nature 
of the transition but the question is still unsettled. 

DNA theories have to account for the fact that, given 
a base pair, the proton forming hydrogen bonds may be 
exchanged between the pair mates whereas proton ex- 
change may not occur in the adjacent base pair. This 
amounts to say that base pair opening, the breathing of 
DNA [24[, is a highly localized phenomenon. Accord- 
ingly the stacking coupling should be weak enough to 
allow for large conformational changes and high flexi- 
bility along the molecule backbone [2^, (2(|. Moreover, 
even when no enzymes or thermo-mechanical forces are 
at work [13, HH , the molecule undergoes large amplitude 
displacements showing the nonlinear nature of the dou- 
ble helix [29N3ll ] . Nonlinearity is intrinsic to DNA and 
should not be attributed to heterogeneity as it persists 
also in artificial homogeneous molecules. Again, fluo- 
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rescence methods have proved that large fluctuational 
openings take place in DNA hairpins while, for double 
stranded DNA, breathing fluctuations induce bubbles of 
a few base pairs which open spontaneously and can be 
monitored in real time [321 ] - Thus dynamical models for 
bubble formation and denaturation should incorporate 
fluctuation effects [U H3| • 

In this regard, a computational method based on the 
imaginary time path integrals [35| has been recently pro- 
posed [36[ to investigate the DNA thermodynamics in 
the temperature range for which denaturation is expected 
to occur. The base pair stretchings with respect to the 
ground state are thought of as time dependent paths with 
the time being an inverse temperature in the spirit of the 
Matsubara's Green functions formalism [37| . The par- 
tition function for a nonlinear PB Hamiltonian is com- 
puted by summing, in the path configuration space, over 
an ensemble of path fluctuations around the ground state 
consistent with the model potential physical constraints. 
The denaturation appears as a rather smooth crossover 
both in homogeneous and heterogeneous DNA fragments 
[38j | and, in the latter, the AT- rich regions drive a multi- 
step transition as a function of temperature. Precisely for 
the same sequences considered in Ref. 38] , a recent anal- 
ysis based on the extended transfer matrix approach 
40] finds consistent results regarding the temperature lo- 
cation of the main peaks in the specific heat profiles. 

Although AT-rich bubbles form at lower temperatures, 
double helix openings extend also well inside the guanine- 
cytosine (GC) domains indicating a role for nonlocal ef- 
fects in shaping multistep denaturation patterns [U IHJ . 
The sequence pattern is particularly relevant in segments 
made of a few tens of base pairs. This short scale, key 
to those transcription starting domains where the genes 
are read, is considered in this study. 

While the path integral method has proved effective in 
dealing with a large number of base pairs displacements, 
modeling the denaturation still presents several unsolved 
questions even for a single molecule of finite length [43j . 
Provided that a molecule state corresponds to a point in 
the configuration space, statistical averages of physical 
quantities should be obtained in principle by summing 
over all states with their Boltzmann weight. In prac- 
tice numerical methods always require a confinement of 
the configuration space 0, HU and, in the path integral 
method, such restriction naturally occurs by selecting the 
ensemble of possible paths for the base pair elongations. 
The size of the ensemble depends however on the model 
potential. 

In this paper, I develop the analysis carried out in Ref. 
[38| and investigate the effects of the backbone stacking 
interactions both in shaping the path configuration space 
for a DNA sequence and in determining the fractions 
of open base pairs versus temperature. The Hamilto- 
nian model is reviewed in Section 2 and the path inte- 
gral method is described in Section 3. The results for 
a short heterogeneous fragment are presented in Section 
4 where, in the denaturation temperature range, I cal- 



culate: a) the fractions of open base pairs, b) the free 
energy derivatives and c) the size of the path ensemble 
by varying the strength of the backbone stacking. The 
confinement of the path amplitudes is also discussed by 
varying the cutoff in the path integration. Some Conclu- 
sions are drawn in Section 5. 



2. Nonlinear Hamiltonian Model 

The nonlinear PB Hamiltonian for a chain of N 
heterogeneous base pairs reads 



N 



^ + V s (y n ,y n - 1 ) + V M (y n ) 

n=l L 1 

K 

V s {y n ,y n -i) = —g(y n ,yn-i){yn ~ y n ~iY 
g(y n ,y n -i) = 1 + pexp[-a(y n + y n -i)] 
V M (y n ) = -D„(exp(-a„y„) - l) 2 , 



(1) 



where y n , the transverse stretching at the n-th site, 
is the degree of freedom for the one dimensional model 
|2j and measures the relative pair mates separation from 
the ground state position. Although Eq. (JTJ does not 
describe the actual geometry of the molecule, it con- 
tains the main contributions which energetically com- 
pete and determine the melting transition, fj, is the re- 
duced mass of the bases which is taken identical both 
for GC- and AT- base pairs [46]. Consistently the back- 
bone harmonic coupling K is site independent: K = p,v 2 
with v being the harmonic phonon frequency. Also the 
non linear (positive) parameters p and a in the stack- 
ing potential Vs{y n ,yn-i) are independent of the type 
of bases at n and n — 1. While the base pair stacking 
interactions generally stabilize the double helix [2^, [26[ , 
I have checked that the homogeneity assumption along 
the molecule backbone has scarce effect on the thermo- 
dynamics in the denaturation range. Instead, the latter 
is essentially determined by the size of K as it will be 
discussed below. Although the form for Vs{y n ,y n -i) in 
Eq. (fTJ may not be unique and other potentials have been 
proposed in the literature pol . |44| , anharmonic stacking 
described by p and a is key to the model as it establishes 
the link with the cooperative character of the bubble for- 
mation. Widely used values are taken in the following 
calculations, p = 1 and a = 0.35A fl9l l47j. 

In fact, when the molecule is closed, t/njj/n-i "C or 1 
for all n, the effective stacking coupling is K(l + p) and 
a large contribution to the stacking interaction comes 
from the overlap of the ir electrons of the base plateaus. 
Whenever either y n > a -1 or y n -\ > a -1 , the corre- 
sponding hydrogen bond breaks and the base moves out 
of the stack thus reducing the electronic overlap. Accord- 
ingly the interaction between neighboring bases along the 
strand weakens and also the next base tends to open. 
Formally, in Eq. (fT]), g(y n ,y n -i) ~ 1 and the effective 
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stacking coupling between neighboring bases drops to K. 
This is the microscopic picture for the interplay between 
anharmonicity and cooperativity which determines the 
formation of a region with open base pairs. When the 
stacking gets weaker, the two strands become more flex- 
ible and their entropy increases. The actual rate of such 
increase approaching denaturation and above will deter- 
mine the character of the transition. 

Hydrogen bonds linking the two bases on complemen- 
tary strands are modeled by the Morse potential Vm(z/?i) 
in Eq. (Q]) 48] which also incorporates the effects of het- 
erogeneous sequences through D n and a n . D n is the dis- 
sociation energy of the pair: for DNA, energies per hydro- 
gen bond are usually taken in the range ~ 15 — 25meV. 
As AT- and GC- base pairs have two and three bonds 
respectively, I take Dat = SOmeV and Dec = 45meV. 
a n is an inverse length defining the spatial scale of the 
potential: transverse stretchings for GC are stiffer than 

for AT base pairs [49|. Accordingly I set a at = 4. 2 A 
and aac — 5A while slightly different values can be 
found in the literature. 

Vm(j/«) is standard for chemical bonds and reproduces 
the main properties of the inter-strand forces: a) it has 
a hard core (y n < 0) corresponding to the repulsion of 
the charged phosphate groups of the backbone; b) it has 
a minimum (y n = 0) for the ground state equilibrium 
distance between the two bases; c) it has a plateau for 
large y n accounting for the fact that the force between the 
bases vanishes at the dissociation energy. Certainly the 
threshold stretching for base pair opening cannot be set 
a priori merely according to the potential parameters, 
it should be rather estimated statistically through the 
mesoscopic model as discussed in Section 4. 

Experiments show that bubbles of only a few base pairs 
are formed and undergo large amplitude thermal fluc- 
tuations which are highly localized I32B . Moreover the 
lifetimes of open and closed states [5fJ are sensitive to 
the solvent as separated strands can recombine at a rate 
which depends on the finite proton concentration in so- 
lution. This recombination event cannot be foreseen by 
the PB Hamiltonian which essentially deals with a DNA 
chain in a infinitely diluted solution. Eq. (TTJ) does not ac- 
count for re-closing as the open strands can go infinitely 
apart with no energy cost due to the plateau in Vm (?/«)• 
Although the model may be improved by changing the 
shape of the potential for base pairs hydrogen bonds , 
a restriction of the configuration space is always required 
to keep the transverse stretchings within a finite range 
40]. This is done in the path integral method as out- 
lined in the next Section. 



by integrating a Boltzmann-like probability distribution 
over the particle path phase space 



Z Q = <j)Dxcxp[-A Q {x}/h] 



(2) 



where Aq{x} is the Euclidean action for the system, 
H is the Planck constant over 2ir, Dx is the integration 
measure and <f indicates that the paths x(r) are closed 
trajectories 52j. r is the imaginary time whose domain 
is r G [0, 0], J3 being the inverse temperature [53j . 

In Ref.[36|, I have proposed a path integral approach 
to the Hamiltonian in Eq. (TTJ introducing the mapping of 
the real space interactions onto the imaginary time scale. 
Accordingly, the transverse stretching y n is represented 
by a one dimensional path x(j~i) obeying the periodicity 
property, xfc) = xfc + (3). In fact there are N + 1 base 
pairs in Eq. (TTJ) . Thus the period j3 can be partitioned in 
N T + 1 points Tj, each of them corresponding to a specific 
base pair n along the DNA strand: 



y n -t x(n), 

n= I ... A ; i= 1 ... N T 

n = ; t Nt+1 ee . 



(3) 



Thus, at any given temperature, the finite size system 
of N+l base pairs, is described by N T + 1 paths taken at a 
specific Ti along the time axis. The presence of an extra 
base pair yo in Eq. (UJ) is remedied by taking periodic 
boundary conditions, yo = yN, which close the finite 
chain into a loop. Such conditions are easily incorporated 
in the path integral formalism as the path is a closed 
trajectory, x(0) = x(/3). Hence a molecule configuration 
is given by N T paths and, in the discrete imaginary time 
lattice, the separation between nearest neighbors base 
pairs is At = (3/N T . Then, the mapping of Eq. ([1]) on 
the time axis also requires: 



y n -i -^x(t') 
t' = Ti - At . 



(4) 



Finally note that the real time derivative y n maps onto 
the imaginary time derivative x(t) according to: 



(5) 



This amounts to replace 



3. Path Integral Method 

In the finite temperature path integral formalism, the 
quantum statistical partition function Zq is viewed as an 
analytical continuation of the quantum mechanical par- 
tition function to the imaginary time and it is calculated 



h^{vpy\ (6) 

which is justified in the classical regime appropriate 
to DNA denaturation. The same replacement is done 
in order to solve the pseudo-Schrodinger equation for a 
Morse potential [54[ which is obtained from Eq. ([TJ in 
the large K limit (and p = 0) [17] . 
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Then, applying Eqs. (|3]) - (J5|) to Eq. (p}, the classical 
partition function for the DNA molecule in one dimension 
reads 

Z c = £z>xexp[-l3Ac{x}] 

Ac{x} = E[f + Vs(x(Ti),x(r')) + V M {x{n)) 

i= 1 

Mf 

x(n) = x Q + [a m cos(w m Ti) + b m sin(cj m 7i) 

ni—l 

uj m = 2m7r//3 




(7) 

where Ac {a;} is the classical Euclidean action. I take 
N T = 100 as the focus is here on short DNA fragments 
[55j . The base pair path x(r^) can be Fourier expanded 
by virtue of the periodic boundary conditions, Mp is de- 
termined numerically and the measure of integration ®x 
is expressed through the coefficients {xq, a m ,b m }. A M is 
the thermal wavelength whose form in general depends 
on the model whether quantum [561 ] or classical. Due to 
Eq. ©, A p = y/lFJjiK. 

The physical features of the path integral approach 
with the associated mapping technique can be synthe- 
sized as follows: 

a) One point (x , a, m ,b m ) in the space of the Fourier 
coefficients selects an ensemble {x(ri),i — 1, N T } repre- 
senting a molecule state with its peculiar base pair dis- 
placements. The state temperature dependence is moni- 
tored by varying j3 as each n e [0,(3]. 

b) Spanning the Fourier coefficients phase space, at 
a fixed (3, one builds a set of possible molecule config- 
urations characterized by different choices of base pair 
stretchings for the same temperature. The higher is the 
number of configurations included in the path integral 
the more reliable is the computation of the DNA ther- 
modynamics. 

c) The truncation of the configuration space mentioned 
in the Introduction occurs in the path integral method 
through the temperature dependent cutoffs A^ and At 
in the latter of Eq. © f57L l58(| . Both A^ and A T are de- 
rived consistently with the requirement that the measure 
of integration normalizes the free particle action and with 
the physics contained in the model potential. This is ex- 
plained in Subsection 4.3 where the possible cutoff model 
dependence of our thermodynamical results is also exam- 
ined. Paths x(ri) ~ represent the equilibrium configu- 
ration for the double helix corresponding to the minimum 
Vm (x( T i))- Larger paths around the minimum should 
be incorporated in the computation however discarding 



negative path amplitudes, smaller than x rn i n ~ — 0.2A 
at high T, which are forbidden by the electrostatic re- 
pulsion between the sugar-phosphate backbones. The 
lower bound confinement for the {x(ri)}, physically due 
to the hard core potential, also ensures that the base pair 
paths are self-avoiding at complementary sites along the 
strands. In fact base pair mates do not overlap. Then 
self- avoidance effects are included in the one dimensional 
model. 

Also positive path amplitudes, larger than x max ~ &A 
at high T, should be discarded posing a threshold on the 
two strands separation. Inclusion of very large positive 
paths does not contribute to the free energy derivatives. 
Instead, an upper bound on the path displacements is 
consistent with strand recombination which may occur 
in the presence of a solvent. Thus the range for the base 
pair stretchings, xfa) € [x m i n ,x max ], is qualitatively set 
according to the shape of Vm^^i)). Inside such range, 
the selection of the paths which indeed contribute to Zc 
is performed by means of a macroscopic constraint, the 
second law of thermodynamics. The strategy is the fol- 
lowing: 

Eq. ([7]) is computed, at an initial temperature Tj, for 
a given path ensemble defined (for any base pair) by the 
number of integration points over the Fourier coefficients, 
N e ff. The latter is temperature dependent. Then, at 
any larger T, the numerical code re-determines the con- 
tribution to Zq and calculates the derivatives of the free 
energy F = — In Zc- 

If, for a given N e ff, the growing entropy constraint 
is not fulfilled then the size of the path ensemble is in- 
creased. The procedure is reiterated until a minimum 
number of paths is found such that the entropy grows ver- 
sus T. This method sets the size of the ensemble whose 
paths satisfy boundary conditions and macroscopic phys- 
ical constraints. N e ff is the T-dependent number of dif- 
ferent trajectories followed by a single base pair stretch- 
ing in the configuration space. As the procedure holds 
for any Tj , the total number of paths contributing to the 
thermodynamics is N T x N e f f . This value sets the over- 
all system size. Numerical convergence has been found 
taking N eff ~ 1.2 x 10 5 at 7/ = 260K while there are 
no significant effects by further increasing the initial size 
of the path ensemble. 



4. Results 
4.1. Fractions of Open Base Pairs 

The path integral method can characterize the melting 
through the computation, via Eq. (|7|) , of the equilibrium 
thermodynamic properties. The latter provide a signa- 
ture of the disruption of the base pair bonds. In fact, 
the order parameter in homogeneous DNA denaturation 
is usually taken as the fraction of bound base pairs ( fb) 
which is proportional to the internal energy of the chain, 
hence —dfb/dT is proportional to the specific heat. Also 
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in heterogeneous sequences the specific heat is an indica- 
tor of the melting as it displays sharp peaks at the tem- 
peratures where various fragments of the sequence open 
[59j | . While the definition of order parameter is less clear 
in this case [l?} , such peaks should be revealed by the be- 
havior of fb ■ Theoretical models should therefore predict 
both macroscopic and microscopic denaturation features 
in the same temperature range although the overlap is 
not expected to be precise for heterogeneous DNA. 

The inciting transition is usually monitored by an in- 
crease in the UV absorption signal around 2600 A. This 
is due to the fact that, when the double helix open and 
the bases move out of the stack, the corresponding elec- 
tronic transitions are less screened. However spectro- 
scopic measurements yield average fractions of open base 
pairs where the averaging is meant over an ensemble of 
molecules. This poses some questions to theorists and 
experimentalists: 

i) There is some intrinsic arbitrariness in the definition 
of an open base pair as one needs to assume a threshold 
for the elongations beyond which the pair dissociates. 

ii) When the UV signal measures that half of the 
base pairs are open, this may indicate either that half 
of the molecules are open and half are closed or that all 
molecules are half-open. Accordingly these methods can- 
not distinguish intermediate states for a single molecule 
configuration which, instead, would be quite interest- 
ing in order to understand the nature of the melting 
transition. New techniques based on quenching of sin- 
gle strands are becoming available to trap intermediate 
states |60| . 

Eq. (J7]) accounts for an ensemble of different configura- 
tions for one molecule. Statistical averages are therefore 
carried out over such ensemble and the mean elongation 
for the i — th base pair reads: 



used to define the fraction of open base pairs / 
The latter is expressed in terms of Eq. © as: 



< x(n) >- 



<j> Dxx(Ti) exp[-/3A c {x}] 



(8) 



where the integral includes those good paths which are 
selected by the method described in the previous Sec- 
tion. Introducing a threshold £ , a base pair is open if: 
< x(ri) > > (. Thus we are able to compute the pro- 
files < x(Ti) > — £ as a function of the site index i for 



different temperatures [61[. C is an arbitrary parameter 
generally depending both on the length and on the se- 
quence of the fragment: for the present model it may 
be reasonably taken [2l|, |4(| in the range ~ [0.5 — 2]A 
thus checking whether, for a given £, there is a signifi- 
cant fraction of base pairs which dissociate close to some 
temperature. Such T values may then be signatures of 
denaturation steps. Too small £ would lead to the wrong 
conclusion that all base pairs are already open at too low 
temperatures. Too large £ would prevent us to follow 
the multistep melting of different regions of the chain at 
different temperatures. 

As the UV signal changes quite abruptly when base 
pairs dissociate, Heaviside step function #(•) is generally 



1 Nr 

-^0(<x(r,)>-C) 



1-/6 



(9) 



Eq. (J9j) is consistent with the definition adopted in 
Monte Carlo simulations [45J and dynamical mesoscopic 
models (62|. 

Alternatively, in the case of a very short experimental 
signal which is not related to < x(ri) >, one may rather 
perform the statistical averages < 8{x(Ti) — C,) > and 
define the fraction of open base pairs /' as 



1 Nr 

f'= ^E(w-o) 



(10) 



While there is no a priori reason to prefer Eq. © or 
Eq. (|10p . I compute here both expressions to verify which 
one is more suitable to capture the multistep denatura- 
tion of an heterogeneous molecule in the framework of 
the path integral formalism. 

I consider a N T =100 fragment with an extended AT- 
substrate in the right side and a slight predominance of 
GC-pairs on the first (from left) 48 sites. The left part 
of the segment has the same sequence as the L48AS frag- 
ment of Ref.[6(| although our model cannot distinguish 
a configuration GC- followed by GC- along the backbone 
from a GC- followed by CG-pair. Experimentally these 
configurations present some differences in the free energy 
stackings [25] ■ The sequence is: 



GC + 6AT + GC + HAT + 8GC + AT + AGC + 
AT + 4GC + AT + 8GC + [49 - 100] AT , 

(11) 

for which bubbles may open in the leftmost and, more 
likely, in the right part of the chain. 

The thermodynamics results are presented in Fig. ^ for 

o —2 

a backbone stiffness K = 60 meVA which corresponds 
to a weak intra-strand coupling. The entropy (Fig. 1(b)) 
is a continuous function of temperature with a slight kink 
at around T — 327K (magnified in the inset) which is 
mirrored by a pronounced peak in the specific heat plot 
(Fig. 1(a)). The specific heat also displays an array of 
minor side peaks at T ~ 275, 300, 380 A". The size of the 
entropy step is ~ 10 meVK , much smaller than the 
value 0.0646meVA~ 1 found in Ref.[l9j], albeit for homo- 
geneous DNA. Even bigger melting entropies have been 
predicted for long chains of heterogeneous DNA using ex- 
tended transfer matrix approach for the PB model but 
with a very large non linear p [63| . Instead, the smooth 
entropy behavior found in the path integral approach is 
more in line with the transfer matrix analysis for a se- 
quence of 100 base pairs of the Joyeux-Buyukdagli model 
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FIG. 1: (Color online) (a) Specific Heat and (b) Entropy 
for the fragment in Eq. (j 1 1 p and stacking coupling K — 
60meVA . (c) Fraction of open base pairs according to 
Eq. ([9|. (d) Fraction of open base pairs according to Eq. (|10|l . 



[44 |64| which also estimates reduced melting entropy 
(with respect to the PB model) at the thermodynamic 
limit. 

Look now at Fig. 1(c) which reports / in Eq. ([9]) for 
three selected values of £: at T ~ 275K, roughly 8% of 
the average base pair stretchings become larger than £ = 
1.5 A, 50% become larger than £ = 1 A and all average 
paths become larger than £ = 0.8 A. At T ~ 295K, 
there is a further enhancement (15% more) in the fraction 
of paths exceeding Q — 1 A while / grows smoothly for 
£ = 1.5 A. At T ~ 327K, all average paths become larger 



than £ = 1 A consistently with the main peak in the 
specific heat. At T ~ 380.fr, / with £ = 1.5A has another 
small increase. The fact that the /— plots display step- 
like enhancements at various T also explains why the 
specific heat peaks are very narrow in temperature. 

Using /' in Eq. (fTDj) . we would not be able to capture 
an evident correspondence between base pair dissocia- 
tions and specific heat peaks: Fig. 1(d) makes clear in 
fact that /' grows steadily for any £. Thus, it seems that 
the averaging method in / is more suitable to account for 
a multistep melting at least for the one-molecule system 
here treated. Our conclusion is in qualitative agreement 
with Refs.[22|,[45| although the investigated sequences dif- 
fer. Instead, in Ref. 21], the expression in Eq. (fT0|) has 
been preferred. Note that, for a given £, /' <C / at any 
T. 

The overall picture emerging from Fig. [T] is that of 
transition driven by the AT-rich portions of the frag- 
ment. Such transition is a smooth crossover (as shown 
by the entropy plot) mainly occurring in the range T ~ 
[275 — 327] K within which all pair displacements become 
larger than £ — 1 A. Certainly, this does not imply that 
the molecule is open at T ~ 327-fT as many more paths, 
mainly associated to the GC-pairs, broaden at higher T 
providing a further entropic gain. In this regard, a recent 
neutron scattering study has pointed out how significant 
clusters of base pairs continue to be closed inside the de- 
naturation regime [65]. In fact the entropy still grows 
above the main denaturation peak although the rate of 
the growth gradually decreases. While we cannot define 
a threshold for the complete strands separation within 
the considered temperature range, we observe that 70% 
of the path stretchings belong to the range [1 — 1.5] A at 
the upper value T = 390K while 30% are larger. 

The calculations for a stacking coupling K = 
100 meVA , are presented in Fig. [5] The entropy plot 
(Fig. 2(b)) is even smoother than in the previous case 
so that the specific heat peaks (Fig. 2(a)) are much less 
sharp. Hardening the backbone stiffness has the main 
effect to shift the melting signatures towards higher tem- 
peratures. Thus, see Fig. 2(c), all average base pairs be- 
come larger than £ = 0.8 A at T ~ 310K, a 35K increase 
with respect to Fig. 1(c). At about the same T there is 
also a 20% increment in the fraction of average paths 
larger than £ — 1 A. For the same value, / shows two 
more step-like features at T ~ and T ~ 385K. 

Again we find an overall correspondence with the spe- 
cific heat peaks locations at T ~ 310!^, T ~ 340i^ 
and T ~ 375K respectively. Also in this case /' grows 
smoothly, see Fig. 1(d), providing no indication of bubble 
formation along the strands. 



4.2. Cooperativity 

Theories for DNA have focussed since longon the co- 
operative character of the melting transition [H, Ho[ . The 
interplay between stacking coupling and cooperativity in 
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FIG. 2: (Color online) (a) Specific Heat and (b) Entropy 
for the fragment in Eq. (j 1 1 p and stacking coupling K — 
WOmeVA . (c) Fraction of open base pairs according to 
Eq. ([9|. (d) Fraction of open base pairs according to Eq. (|10|l . 



the model of Eq. (JTJ has been pointed out in Section 2. 
In the path integral method the degree of cooperativity 
is measured by the parameter N T x N e f t which yields the 
total number of trajectories contributing to Eq. ([7]) at a 
given T . This is set by the numerical code as described in 
Section 3. I emphasize that N T x N e f / represents (for any 
T) the minimum number of paths for which the second 
law of thermodynamics is fulfilled. Once such minimum 
value is determined there is no need to further increase 
the size of the paths ensemble. In Fig. [31 N T x N e f f is 
plotted versus T for the sequence in Eq. as a function 



FIG. 3: (Color online) Total number of paths used in the 
computation of the thermodynamics of the DNA sequence 
for five values of the stacking parameter K in units meVA 



of the stiffness K. 

While the computation starts with 1.2 x 10 7 paths for 
any K at the lower bound of the T range, the renormal- 
ization of N T x N e f f along the T-axis remarkably depends 
on K. At T — 390X, about 3 x 10 7 paths are included in 

o —2 

Zc for the case K — 60meVA whereas such number 

o —2 

drops by a factor two in the case K = 100 meVA . For 
the former case, also note the step-like enhancement at 
T ~ 327 K consistent with the main peak in Fig. 1(a). 
An enhanced stiffness reduces the flexibility of the strand 
and inhibits the bubble formation thus shifting the oc- 
currence of the melting steps at higher T. Here we see 
how the intra-strand stacking is key to determine the 
inter-strands opening probabilities. The lowest K values 
here assumed are in the range of those usually taken in 
nonlinear PB Hamiltonian models [l9j |. 



4.3. Path Amplitudes Cutoffs 

The cutoffs Ay and At in Eq. permit to build an 
ensemble of paths whose amplitude is temperature de- 
pendent. To derive analytically the form of the T— de- 
pendence I observe that the integration measure normal- 
izes the free particle action: 



<j)T)x(T)cxp — J (It~x(t) 



(12) 



By inserting the Fourier path expansion in Eq. ([?)), the 
l.h.s. of Eq. (fT2j) transforms into a product of Gaufiian 
integrals which yield the mathematical criteria to set A^, 
and At through the conditions: 
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—= / ds m exp(-s^) = 1 
Vf Jo 



2 _ 



\2 

l-i 



(13) 



with U being a dimensionless cutoff which can be set 
numerically by taking, for instance, a series expansion 
for the Gaufiian integral 66]. From Eqs. ([T3| . it follows 
that: 



A° 



A T = 



<7A 



r 3/2' 



(14) 



T be- 



hence, the path amplitude cutoffs display a oc 
havior. 

There is however some freedom in the choice of U, that 
is one may tune the cutoff within a range which essen- 
tially fulfills both Eq. (fT3"|) and the physical requirements 
of the problem in Eq. (J7J emphasized in the previous Sec- 
tion. In this way we have a mathematically consistent ap- 
proach to explore the hydrogen bonds potential plateau 
and include a set of base pair paths suitable to the spe- 
cific temperature of the system. All the results presented 
so far have been obtained taking U — 17 that implies 
A389.R" = 4.04A for the first Fourier coefficient. Accord- 
ingly the path displacements are not larger than ~ 5 A in- 
side the investigated temperature range. Fig. EJa) shows 
the total number of paths larger than 1, 2, 3, 4 A respec- 
tively, used in the computation of Fig. Q] Note that a 
significant number of paths, ~ 10 5 , larger than 4 A but 
smaller than 5 A contributes to Zq at T ~ 350A". 

How do the thermodynamical properties depend on the 
choice of U? Taking U = 22, which means A^ggfc — 
5.23A, I incorporate even broader base pair displace- 
ments as made evident by Fig. Hfb) : at T ~ 350 AT there 
are now about 10 5 paths, larger than 5 A but smaller than 
6 A. Specific heat and entropy are plotted in Fig. @|c) 
and Fig. EtA) respectively. A remodulation in the specific 
heat peaks occurs with respect to Fig.QJa): three sizeable 
side peaks (rather than two) show up in the specific heat 
at temperatures lower than the main peak which is now 
located at T — 342AT. Some portions of the fragment 
tend to open before the main transition temperature is 
attained while the latter is shifted slightly upwards in 
comparison with the case of Fig. QJa). As a main result, 
the entropy jump at T = S42K is again small and of 
order 10~ 4 meVK~ 1 , see inset in Fig.^Jd). Likewise tak- 
ing U — 100, A 3 89K = 23.79A, thus confirming that the 
smooth character of the denaturation transition is not an 
artifact of the truncation procedure in the path config- 
uration space. In the path integral method, the cutoff 
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FIG. 4: (Color online) (a) Number of paths larger than 
1,2, 3, 4 A taking in Eq. the same cutoff, 4.04 A at T = 
389^, as in Fig. [T]and [2] (b) Number of paths larger than 
1,2,3,4, 5 A assuming a larger cutoff, 5.23 A at T = 3897f, 
in Eq. (ff)l. (c) Specific Heat and (d) Entropy calculated with 
the At value given in (b). K in units meVA 



dependence appears altogether not so crucial as the en- 
semble size N T x N e ff is large even for short sequences. 



5. Conclusion 

I have applied the path integral formalism to an hetero- 
geneous DNA sequence modeled through the nonlinear 
Peyrard-Bishop Hamiltonian which captures the main 
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interactions of the double strands configuration. The 
method allows us to compute the thermodynamic prop- 
erties incorporating all base pairs thermal fluctuations 
that are particularly relevant in a short fragment. The 
energetic gain associated to the (bounded) double strands 
molecule competes with the entropic gain due to the large 
number of configurations available once the two strands 
begin to dissociate. The computation refers to a temper- 
ature window in which denaturation is expected to occur 
and includes an high number of possible trajectories for 
the base pair stretchings. Such number is set at any 
temperature on the base of physical criteria included in 
the numerical code. First, the path configuration space 
is truncated consistently with the requirements of the 
hydrogen bonds model potential. Second, the allowed 
base pair paths are selected by simply requiring that the 
entropy should have a positive temperature derivative. 
Such requirement does not introduce any constraint on 
the shape of the entropy itself which always appears as 
a continuous function of T. Thus the denaturation man- 
ifests as a smooth crossover taking place in multisteps 
as the average stretchings of the adenine-thymine base 
pairs are larger than those of the guanine- cy to sine base 
pairs. 

The work has focused on the quantitative effects of the 
backbone stiffness and, in particular, it has been shown 
that such parameters determines the size of the path 
ensemble included in the computation of the partition 
function. This points to a relevant role of the stacking 
coupling on the degree of cooperativity of the molecule. 
A weak stacking induces flexibility of the double strand 
structure and allows an high number of paths to partici- 



pate to the multistep melting. The importance of anhar- 
monicity in the intra-strand potential with regard to the 
melting transition is controversial and has been widely 
discussed in the literature: we recognize that the anhar- 
monic stacking drives a cooperative behavior of the base 
pairs in the PB Hamiltonian model but the effects of the 
coupling on the denaturation features are due to the har- 
monic parameter K rather than to the anharmonic force 
constants. It is the value of K which sets the location of 
the melting steps along the T— axis and determines the 
fractions of open base pairs. However varying the stack- 
ing coupling does not change the character of the transi- 
tion which remains smooth as in homogeneous fragments. 
In particular, I find that the entropy jump associated to 
the main melting transition is very small due to fluctua- 
tional effects which are strong in a short chain, whereas 
the fractions of open base pairs display step- like enhance- 
ments in remarkable correspondence with the peaks of 
the specific heat. This investigation also sheds light on 
the averaging procedure to estimate such fractions for a 
molecule existing in many different configurations. 

Among the main advantages of the path integral 
method there is certainly the capability to include a 
large-size ensemble of base pair fluctuations within a self- 
contained computational time. While this study refers 
to a one-dimensional system and cannot account for the 
helicoidal structure of the molecule, the method seems 
promising also for extensions to the three dimensional 
space. This would also permit to fully incorporate ex- 
cluded volume effects which have been here only partially 
considered. 
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